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ABSTRACT 

Using analytic arguments and numerical simulations, we examine whether 
chondrule formation and the FU Orionis phenomenon can be caused by the 
burst-like onset of gravitational instabilities (Gls) in dead zones. At least two 
scenarios for bursting dead zones can work, in principle. If the disk is on the 
verge of fragmention, Gl activation near r ~ 4 to 5 AU can produce chondrule- 
forming shocks, at least under extreme conditions. Mass fluxes are also high 
enough during the onset of Gls to suggest that the outburst is related to an FU 
Orionis phenomenon. This situation is demonstrated by numerical simulations. 
In contrast, as supported by analytic arguments, if the burst takes place close to 
r ~ 1 AU, then even low pitch angle spiral waves can create chondrule-producing 
shocks and outbursts. We also study the stability of the massive disks in our 
simulations against fragmentation and find that although disk evolution is sen- 
sitive to changes in opacity, the disks wc study do not fragment, even at high 
resolution and even for extreme assumptions. 

Subject headings: accretion, accretion disks - hydrodynamics - instabilities - 
planetary systems: protoplanetary disks - solar system: formation 

1. Introduction 

1.1. Gls, MRI, FU Orionis Outbursts, and Chondrules 

Gravitational instabilities can activate in a disk when the Toomre (1964) parameter 
Q — Csk/ttGY: < 1.7 for a thick disk (see Durisen et al. 2007), where Cg is the sound speed, 
E is the surface density, and k is the epicyclic frequency. As indicated by Q, a disk is 
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unstable against GIs when it is cold and/or massive. The resulting spiral waves driven by 
self-gravity efficiently transfer angular momentum outward and mass inward (e.g., Lynden- 
Bell & Kalnajs 1972; Durisen et al. 1986). Another mechanism that can efficiently transfer 
angular momentum outward is the magnetorotational instability (MRl; see Balbus & Hawley 
1991; Desch 2004). In contrast to GIs, the MRI only requires a weak magnetic field coupled 
to the gas. These mechanisms, either separately or in combination, are likely to be the 
principal way T Tauri stars accrete gas from a disk (Hartmann et al. 2006). 

In order for the MRI to occur, ionized species must be present in the gas phase. Thermal 
ionization of alkalis occurs wherever T > 1000 K, but depletion of ions by dust grains may 
move the temperature threshold closer to T ~ 1700 K (Desch 1999; Sano et al. 2000), where 
the dust sublimates completely. Elsewhere, the ionization must be driven by a nonthermal 
source, e.g., energetic particles (EPs). For this discussion, EPs refers to any particles that 
could ionize the gas, e.g.. X-rays. Gammie (1996) proposed that disks may have active and 
inactive MRI layers due to attenuation of EPs by the gas. In the inner regions of a disk (at 
radii ~ few AU) where the column densities are large, MRI may only be active in a thin layer, 
resulting in layered accretion. As one moves outward, column densities drop, and the entire 
disk can become MRI active. The region where the MRI is mostly absent is called the dead 
zone. EPs are attenuated by a surface density of only about 100 g cm~^ (Stepinski 1992), 
and so even a minimum mass solar nebula (MMSN) will likely exhibit layered accretion 
(Desch 2004). Even if mass accretion is only reduced and not altogether halted as a result 
of Reynolds stresses (Fleming & Stone 2003; Oishi et al. 2007), mass may still pile up in the 
dead zone. If enough mass accumulates, then even for an otherwise low-mass disk, GIs can 
activate. 

The FU Orionis phenomenon is characterized by a rapid (1-lOs yr) increase in optical 
brightness of a young T Tauri object, typically by 5 magnitudes, and is driven by sudden 
mass accretion of the order 10"'^ yr~^ from the inner disk onto the star (Hartmann & 
Kenyon 1996). Because FU Ori objects appear to have decay timescales of about 100 yr, 
the entire mass of a MMSN (~ 0.01 Mq) can be accreted onto the star. To date, the best 
explanation for the optical outburst is a thermal instability (e.g.. Bell & Lin 1994; Kley & Lin 
1999; see discussions in Hartmann & Kenyon 1996, Green et al. 2006; and Zhu et al. 2007). 
Armitage et al. (2001) suggested that GIs in a bursting dead zone might be able to trigger 
an FU Ori outburst by rapidly increasing the accretion into the inner disk (< 0.1 AU) and 
initiating an MRI through thermal ionization. Hartmann (2007, in private communication) 
and Zhu et al. (2007) also suggest that the heating due to gravitational torques might drive 
an MRI for r > 0.1 AU, which would then feed mass inside 0.1 AU until a thermal instability 
sets in. The FU Ori phenomenon may be a result of a cascade of instabilities, starting with a 
burst of GI activity in a dead zone, followed by accretion due to an MRI, followed finally by 
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a thermal instability (cf Kley & Lin 1999). Indeed, recent observations of FU Ori indicate 
that very large mass fluxes are present out to at least r ~ 0.5 AU (Zhu et al. 2007). 

Although details are still being debated, most chondrules formed in the first 1 to 3 Myr of 
the Solar Nebula's evolution (Bizzarro et al. 2004; Russell et al. 2005). Chondrule precursors 
were flash melted from solidus to hquidus, where high temperatures T ~ 1700 K were 
experienced by the precursors for a few minutes. The melts then cooled over hours, with the 
actual cooling time depending on chondrule type. Chondrule collisional histories and isotopic 
fractionation data, chondrule-matrix complementarity, fine-grained rim accumulation, and 
petrological and parent body location arguments (Krot et al. 2005) suggest that chondrules 
formed in the Solar Nebula (Wood 1963) in strong, localized, repeatable heating events. The 
shock wave model for chondrule formation can accommodate these observational constraints 
and reproduce heating and cooling rates required to form chondrule textures (lida et al. 2001; 
Desch & Connolly 2002; Cielsa & Hood 2002; Miura & Nakamoto 2006). One plausible source 
of chondrule-producing shocks is a global spiral wave (Wood 1996). Marker & Desch (2002) 
suggest that spiral waves could also explain thermal processing at distances as large as 10 
AU, which may be necessary to explain observations of comets (Wooden et al. 2005) and 
recent Stardust results (e.g., McKeegan et al. 2006). It has been suggested that bursts of 
GIs may be able to produce the required shock strengths (Boss & Durisen 2005) and provide 
a source of turbulence and mixing (Boss 2004b; Boley et al. 2005). Global spiral shocks 
are appealing because they fit many of the constraints above. They may be repeatable, 
depending on the formation mechanism for the spiral waves; they are global, but produce 
fairly local heating; they can form chondrules in the disk; and they can work in the inner 
disk as well as the outer disk. 



1.2. Fragmentation 

Knowing under what conditions protoplanetary disks can fragment is crucial to under- 
standing disk evolution inasmuch as a fragmented disk may produce gravitationally bound 
clumps. This has become known as the disk instability hypothesis for the formation of gas 

giant planets (Kuiper 1951; Cameron 1978; Boss 1997, 1998). The strength of GIs is regu- 
lated by the coohng rate in disks (Tomley et al. 1991, 1994; Pickett et al. 1998, 2000, 2003), 
and if the cooling rate is high enough in a low-Q disk, a disk can fragment (Gammie 2001). 
Gammie quantified that a disk will fragment when tcooi^ ^ 3 for a disk with a F = 2, where 
F is the two-dimensional adiabatic index, such that J pdz = P ~ E'", where p is the gas 
pressure and z is the vertical direction in the disk. Here, icooi is the local cooling time and 
Q is the angular speed of the gas. This criterion was approximately confirmed in 3D disk 
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simulations by Rice et al. (2003) and Mejia et al. (2005). Rice et al. (2005) showed through 
3D disk simulations that this fragmentation criterion depends on the 3D adiabatic index 
Fi and, for Fi = 7 = 5/3 or 7/5, the fragmentation limit occurs when icooi^ ^ 6 or 12, 
respectively. These results show that a change by a factor of about 1.2 in 7 has a factor of 
two effect on the critical cooling time. In addition, these results indicate that the cooling 
time must be roughly equal to the dynamical time of the gas for the disk to be unstable 
against fragmentation when 7 = 5/3. 

Do such prodigious cooling rates occur in disks when realistic opacities are used with 
self-consistent radiation physics? This question is heavily debated in the literature (e.g.. 
Nelson et al. 2000; Boss 2001, 2004a, 2005; Rafikov 2005, 2007; Boley et al. 2006, 2007b; 
Mayer et al. 2007; Stamatellos & Whitworth 2008). The simulations to date use a wide 
variety of numerical methods, including very different approximations for radiation physics, 
and only two groups have published results of radiative transfer tests appropriate for disks 
(Boley et al. 2007b; Stamatellos & Whitworth 2008). 

Nelson et al. (2000) used 2D SPH simulations with radiation physics to study proto- 
planetary disk evolution. Because their simulations were evolved in 2D, they assumed that 
the disk at any given moment was in vertical hydrostatic equilibrium. Using a polytropic 
vertical density structure and Pollack et al. (1994) opacities, they cooled each particle ac- 
cording to an appropriate effective temperature. In their simulations, the cooling rates are 
too low for fragmentation. In contrast. Boss (2001, 2005) employed radiative diffusion in his 
3D grid-based code; fragmentation occurs in his simulated disks. Besides the difference in 
dimensionality of the simulations. Boss assumed a fixed temperature structure for Rosseland 
mean optical depths less than 10, as measured along the radial coordinate (cf recent flux- 
limited simulations by Boss 2008). Boss (2002) found that the fragmentation in his disks is 
insensitive to the metallicity of the gas and attributed this independence to fast cooling by 
convection (Boss 2004a). However, it must be noted that Nelson et al. (2000) assumed a ver- 
tically polytropic density structure. Because the specific entropy s ~ In where p — Kp^ 
is the polytropic equation of state and p is the gas density, the Nelson et al. approximation 
effectively assumes efficient convection. 

Except for extremely massive and extended disks (Stamatellos et al. 2007), recent simu- 
lations with radiative physics by Cai et al. (2006), Boley et al. (2006, 2007b), and Stamatellos 
& Whitworth (2008) find long cooling rates and no efficient cooling by convection. Cai et 
al. also show that the strength of Gls is dependent on the metallicity. Furthermore, Boley 
& Durisen (2006) suggest that shock bores, which can cause a rapid vertical expansion in 
the post-shock region of spiral shocks, could be misidentifled as "convective" flows. In both 
contrast and support of these studies, Mayer et al. (2007) use SPH simulations with 3D flux- 
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limited diffusion and find that fragmentation only occurs once the mean molecular weight of 
the gas fi > 2.7. However, the simulations presented by Meji'a (2004), Cai et al. (2006), and 
Boley et al. (2006) unintentionally were evolved with fi ^ 2.7 due to an error in the inclusion 
of He in the opacity tables, and their disks do not fragment. The issue of fragmentation in 
radiatively cooled disks thus remains unsettled. 



1.3. Current Study 

For this study, we adopt the hypothesis, which we refer to as the unified theory, that 
bursts of GI activity in dead zones drive the FU Ori phenomenon and produce chondrule- 
forming shocks. This hypothesis is an amalgamation of ideas presented in Wood (1996, 2005), 
Gammie (1996), Armitage et al. (2001), and Boss & Durisen (2005). In order to investigate 
this scenario, we designed a numerical experiment to evolve a massive, highly unstable disk 
with an initial radial extent between 2 and 10 AU. Commensurately, we investigate the 
stability of these massive, gravitationally unstable disks against fragmentation to assess the 
feasibility of the disk instability hypothesis for gas giant formation. 



2. Expectations 

The shock speed Ui for a fluid element entering a global, logarithmic spiral wave with 
pitch angle i is approximately described by 

.,.30kms-(-j (-) |l-y |sin, (1) 

where is the corotation radius of the spiral pattern and where we have assumed that the 
gas azimuthal motion is Keplerian and Vr is negligible. When ^ r, the shock speed limits 
to the Keplerian speed times sini; whether spiral waves with large Vp produce chondrules 
is mainly dependent on the pitch angle of the spiral wave. In contrast, when r ^ Vp, the 
speed increases as r, and even shallow pitch angles can produce chondrules. For example, 
consider the MMSN, where the midplane p(r) = 1.4 x 10"^ (r/1 AU)~"^^^^ g cm~^ (Hayashi et 
al. 1985). Figure 1 shows the Ui-p plane with heavy, solid curves indicating shock speeds for 
Tp = 1, 2.5, 5, and oo AU and with i = 10 and 30°. The colored regions on the plot highlight 
where chondrule formation is expected in a MMSN (yellow) and a lOxMMSM (orange). 
The chondrule-forming curves are based on the results of Desch & Connolly (2002), which 
we summarize by mi ~ — (11 + 21ogp (g cm~'^)) ih 0.5 km s^^. The boundaries set by these 
curves are meant to be illustrative, not definitive. Spiral shocks along shallow pitch angle 
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spiral waves {i ~ 10°) are mostly if not entirely out of the chondrule-forming range between 
r = 1 and 5 AU for shocks inside corotation. However, spirals with corotation in the inner 
disk can produce chondrule-forming shocks for a wide range of radii, depending on the actual 
Tp. If the spiral waves have very open pitch angles {i > 30°), then spiral waves with corotation 
near r = 5 AU can produce chondrule-forming shocks near the asteroid belt and at cometary 
distances. The mass of the disk does not change these general behaviors. A major caveat 
for this simple-minded approach is that the fluid elements' motions may be poorly described 
by equation (1) if vertical and radial excursions induced by shock bores (Boley & Durisen 
2006) cannot be neglected. 

The possibility of producing chondrules in the asteroid belt and at comet distances 
during the same burst is attractive. Moreover, we would also like to use these simulations 
to study other phenomena such as radiation transport, convection, and disk fragmentation 
in the 2-10 AU region. For these reasons, we tried to design a numerical experiment to 
investigate the connection between Gls. FU Orionis events, and chondrule formation that 
is biased toward producing strong shocks with a corotation near 5 AU and with open pitch 
angles. 



3. Methodology 

3.1. Initial Conditions 

If chondrulc-producing shocks cannot be created in simulations biased in their favor, 
then it would present a serious problem for Gls as the source of chondrule processing. The 
reader should keep in mind that the model, at this stage, is not necessarily meant to be 
representative of a typical disk, although the massive disk we present may be plausible 
for early Class I objects (e.g., L1551; Osorio et al. 2003). To create the initial conditions 
(ICs), consider a disk that is vertically polytropic, is axisymmetric, has a constant 7, and is 
Keplerian {k — D,). Also assume that self-gravity is negligible. With these assumptions, the 
vertical disk structure can be calculated by using the following equations: 

piz) = po{l-zyh^f^'^-'\ (2) 

^ - 0^(731)^0 , and (3) 

,^n^7-lV^' r[(37-l)/(27-2)] y^^"^'^ 

= l^H^J r[7/(7-i)] j ' 
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where K for p = Kp^ is the poljd;ropic coefficient at any given r, G is the gravitational 
constant, VL is the Keplerian angular speed, and h is the disk scale height. Using equations 
(2) through (4), one can show that 



To select a model, one needs to specify power laws for S and Q, S and or K and Q. 
For calculating Q, the midplane sound speed Cg is used. Note that when Q = constant, 
there exists a lower limit for g, S ~ r~'^, where the radial entropy gradient remains positive. 
Because the specific entropy s ~ InK, equation (5) can be solved to show that ds/dr ~ 
(37 — 2q)/r for constant Q. The critical surface density power law is then Qc = 5/2 for 
7 = 5/3 and qc = 21/10 for 7 = 7/5. When q > qc, radial thermal convection may set in for 
a constant Q disk. However, the radial entropy gradient is not the only stability criterion 
for radial thermal convection in disks (Klahr & Bodenheimer 2003). 

Equations (2) through (5) are used to generate the ICs for this study. The disk is 
massive, 0.09 Mq, with S ~ r^^, initially between about 2 and 10 AU in radius and with 
a Q ~ 2 everywhere. Although the surface density power law is steep, it is consistent with 
the Nice-MMSN (Desch 2007)0 Because the analytic model ignores self-gravity, assumes 
a constant 7 (see §3.2 below), and has sharp disk truncation edges, the model must be 
relaxed to a new equilibrium structure in the hydro code. Once the ICs are generated, 
they are loaded into our code and evolved without cooling at very low azimuthal resolution, 
essentially two-dimensionally. The radial and vertical momenta are damped to allow the disk 
to relax calmly to an equilibrium configuration, which takes about ten orbits at 10 AU. To 
create a mass concentration in the disk, mass is then added to the disk linearly in time with 
a Gaussian profile. The specific internal energy is held constant. The peak of the Gaussian 
is centered on 5 AU. The radial FWHM of the Gaussian is chosen to be 3 AU, which is 
roughly the size of the most unstable radial wavelength ~ 2ttCs/Qk, ^ 27ih/Q ^ 0.27rr 
(Binney & Tremaine 1987). Mass is added until nonaxisymmetry is visible in midplane 
density images, which takes approximately 190 yr, i.e., 6 orp of evolution (1 orp is defined 
as 1 outer rotation period at 10 AU for this model, which is about 33 yr). The total mass 
added is about 0.08 Mq, and so if one were to imagine a corresponding accretion rate, it 
would be about 4 x 10"'* Mq yr~^. It is probably unphysical to add so much mass so quickly 



^This model for the MMSN is derived in a way similar to Weidenschilling's (1977) approach; but, instead 
of using the current positions of Solar System planets, Desch uses initial positions based on the Nice model 
(e.g., Levison et al. 2007). 
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to the disk without increasing the temperature. However, we remind the reader that the 
study is a numerical experiment. Figures 2 and 3 illustrate the density distributions of the 
disk before and after the mass accretion. 

The mass buildup increases the total mass of the simulated disk to ~ 0.17 and it 
decreases Q below unity in a narrow region around 5 AU just before strong GI activation. 

Because the ring is grown over 6 orp, the instability is not obviously overshot. The mass- 
weighted average Qav over the FWHM of the Gaussian centered at 5 AU is approximately 
unity when the ring becomes unstable. This implies that Q must approach unity for a spatial 
range of at least the most unstable radial wavelength before GIs will activate in a disk. 

3.2. CHYMERA 

The disk is evolved with the CHYMERA code (Boley 2007), which employs the Boley 
et al. (2007b) radiation physics algorithm. CHYMERA is an Eulerian code that solves 
the hydrodynamics equations on a fixed, cylindrical grid. The code includes self-gravity, 
a ray+fiux-limited diffusion hybrid scheme, and a variable first adiabatic index Fi. This 
code has been extensively tested. In response to criticisms of Boss's (2001, 2002, 2004a, 
2005) simulations by Boley et al. (2007a), Boss (2007) questions the accuracy of the Indiana 
University code, which is the predecessor of CHYMERA. Boss raises concerns that we address 
here. 

Potential Solver: CHYMERA uses a direct Poisson solver to calculate the potential due 
to mass on the grid. To provide the boundary conditions for the direct solver, a spherical 
harmonics expansion for Z = |m| < 10 is used to calculate the potential on the grid boundary 
(e.g., Pickett et al. 2003; Boley et al. 2007b) . Boss argues that the boundary solver used 
in CHYMERA may be too low-order to capture clump formation. Tests by Pickett et 
al. (2003) and Mcjia et al. (2005) and the results from the Wengen Test 4 code comparison 
(WT4; Mayer et al. 2008, in prep.) demonstrate that this is not the case. During the 
initial analyses of WT4, the potential solver was rigorously tested in an attempt to explain 
some deviations noticed between CHYMERA and other codes. These deviations were finally 
attributed to differences in the initial conditions, and new simulations with consistent ICs 
show very good agreement between codes. To test the potential solver in CHYMERA, a 
high-order Bessel function expansion boundary solver (Cohl & Tohline 1999) was employed 
with Legendre half-integer polynomials expanded out to m = 30, which showed convergence 
to better than 10~^ when compared with m = 40. The spiral structure and clump formation 
in the simulation with this high-m boundary solver were indistinguishable from the normal 
lower-order spherical harmonics expansion. The Wengen test is very demanding, with strong 
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nonaxisymmetric structure and a highly flattened disk, and yet the lower-order spherical 
harmonics expansion for Z = \m\ < 10 used in CHYMERA (e.g., Boley et al. 2007b) captures 
clump formation. 

Resolution & Numerical Heating: The resolution for the high-resolution simulations 
presented in this paper are at higher resolution than the base simulations presented for 

Wengcn Test 4. As will be demonstrated here, CHYMERA shows fragmentation when it 
is expected based on experimentally determined fragmentation criteria (Gammie 2001; Rice 
et al. 2005). The numerical heating reported by Boley et al. (2006) affected a small region 
of their inner disk, and is understood to be a resolution effect. The simulations presented 
here are at high enough resolution, as determined by the Boley et al. (2007b) tests, that this 
effect should be negligible if present at all. 

Radiative Transfer Boundary Conditions: CHYMERA has been shown to relax to ana- 
lytic solutions for the flux through an atmosphere as well as the temperature profile for high, 
moderate, and low optical depths. The code allows convection when it should, and does not 
when it should not (Boley et al. 2007b). 

For the treatment of H2, we use a frozen ortho/para ratio of 3:1 in this study for 
reasons presented in Boley et al. (2007a). The grain size distribution is assumed to follow 
the ISM power law, dn ~ a~^-^da (D'Alessio et al. 2001), where the maximum grain size 
Omax is chosen to be 1 mm to account for grain growth (D'Alessio et al. 2006). To transition 
smoothly between the low and high optical depth regimes, a weighted combination between 
Rosseland and Planck mean opacities is used Kw — {krTr + Kp / Tr) / {tr -\- ^/tr), where 
tr is the vertical midplane Rosseland optical depth. Different methods for interpolating kr 
and Kp along the ray are being explored, but they have not been implemented here. The 
particles are assumed to be well-mixed with the gas. The effect of dust settling is heuristically 
considered by evolving the disk in several ways: standard opacity and 10~^, 10~^, and 10^^ of 
the standard opacity, referred to as the 10^^ 10"'^ and lO"'' n simulations. This choice 
in opacities varies the vertical midplane Rosseland mean optical depths from about r = 10^ 
to unity, and spans the limits where advection of photons becomes important {rv/c > 1; 
Krumholz et al. 2007) and where radiative cooling becomes most efficient (r ~ 1). No 
external radiation is assumed to be shining onto the disk, except for a 3 K background 
temperature. Although such rapid dust settling and a naked disk are likely to be unrealistic, 
the assumptions are made to bias the disk toward strong shocks inasmuch as irradiation and 
opacity affects Gl strength (Cai et al. 2006, 2008). For a base comparison, an adiabatic 
simulation, i.e., where shock heating is allowed but not cooling, is evolved for about 33 yr. 

The simulations are evolved at two resolutions: (r, 0, z) = 256 x 512 x 64 and 512 x 
1024 X 128 cells, respectively, above the midplane. The lower resolution simulations (512 
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sim) have a grid spacing of 0.05 AU per cell in r and z and the higher 0.025 AU (1024 sim). 
For both resolutions, rA0 ^ Ar at about r = 4 AU. These parameters are summarized in 
Table 1. For each disk, 1000 fluid tracer elements are randomly distributed in r between 
about 3 and 7 AU, in over 27r, and in z roughly within the scale height of the disk. 

Due to inefficient cooling in the standard simulations, high enough temperatures in 
the disk are reached after 80 and 55 yr for the 512 and 1024 simulations, respectively, to 
create radiative timescales that arc too short to resolve; the simulations are stopped. As 
described below, the reduced opacity simulations are able to cool much more efficiently, 
and the disk temperatures remain manageable with the time-explicit algorithms used in 
CHYMERA. However, the reduced opacity simulations are stopped after about 110 yr due 
to the development of a strong, one-arm spiral, which is treated incorrectly with our fixed- 
star assumption. 



3.3. Fluid Element Tracer 



The hydrodynamics scheme in CHYMERA is Eulerian, and does not give direct infor- 
mation on the histories of tracer fluid elements. In order to derive detailed and statistical 
shock information and to capture the complex gas motions in unstable disk simulations, we 
have combined a tri-Akima spline interpolation algorithm with a Runge-Kutta integrator 
(e.g.. Press 1986) for tracing a large sample of fluid elements during the simulation. Ve- 
locities and thermodynamic quantities, e.g., temperature and density, are interpolated once 
each time step. Because the hydrodynamics code explicitly solves the equation of motion for 
the gas, the algorithm only needs time and velocity information to integrate the positions 
of the fluid elements. The Runge-Kutta integrator was intended to be fourth-order accu- 
rate; unfortunately, an implementation error when combining the integration scheme into 
CHYMERA resulted in a first-order accurate scheme in time. Because tracing fiuid element 
histories is an important point of this paper, we reran a stretch of one of the simulations 
with a second-order fiuid element tracing routine. After about one orbit at r ~ 5 AU, the 
mean fractional differences, {x\ — x-i)lx2i between the first and second-order schemes are 
—3 X 10~^, —2 X 10~^, and —10^^ for r, 0, and z. respectively. The sample standard de- 
viations are 4 x 10^^, 10^^, and 10~^ for r, 0, and z^ respectively. The mean differences 
are small, and show that there are no obvious systematic effects. The sample deviations 
indicate that the variations between the the two schemes arc also marginal for the radial 
and azimuthal directions. There is considerably larger scatter in the vertical direction, where 
spiral shocks can create complicated vertical fiows (Boley & Durisen 2006). Because the fiuid 
elements are used to gather shock information, the scatter in the vertical positions between 



-li- 



the two schemes is not a major concern for this study. The effect is similar to redistributing 
the fluid elements vertically, keeping the same r and 0, every few orbits. 

An Akima spline is similar to a natural cubic spline, but typically yields better results for 
curves with sudden changes (see Akima 1970), as are expected in shock profiles. Although the 
spline fits a curve to a one-dimensional set of data points, the interpolation can be extended 
to data in three dimensions. First, consider a cubic volume with data at the vertices. A 
value anywhere in the volume can be approximated by seven linear interpolations: four to 
calculate values between each vertex in a particular direction, two to calculate the values 
along the projections of the desired point onto the interpolated lines, and a final interpolation 
through the point of interest. Extending this from a simple tri-linear fit to a tri- Akima spline 
is relatively straightforward. Instead of using the eight nearest points that enclose a volume, 
one uses the 125 closest points, with five data points used for every Akima spline fit. The 
central data point is the closest data value to the point of interest. We use the GNU Scientific 
Library Akima spline algorithm for performing fits. 



4. Results 
4.1. Structure 

As these disks become gravitationally unstable, low-order spiral waves develop. These 
waves drive a sudden increase in mass flux throughout the disk and drive strong spiral shocks. 
Following Mejia et al. (2005), we refer to this phase of disk evolution as a burst of GI activity. 
Surface density plots for the 512 simulations are shown in Figures 4 and 5 at t ~ 33 and 
77 yr, respectively. The images at t ~ 33 yr indicate that except for the 10""^ k, simulation, 
the burst has a well-defined, three- arm spiral and that the spiral arms become denser as the 
opacity is lowered. By 77 yr, the disk structures are noticeably different. Each disk has a 
visually distinct two-arm spiral, again, except for the 10~^ k simulation. For comparison 
with Figures 4 and 5, the 1024 simulations are shown in Figures 6 and 7 at i ~ 33 and 78 
yr, respectively. The higher resolution simulations have additional fine structure, but they 
do not develop clumps or dense knots. In fact, some of the knotty structure that develops 
in the 512 10^^ k simulation is absent at higher resolution. The lower opacity simulations 
appear to have stronger amplitudes, with denser spiral arms and stronger mid-order Fourier 
components (see below). As will be discussed, the 10~^ k, simulation is on the verge of 
fragmentation. 

The visual differences between the simulations can be quantified by a Fourier mode 
analysis. Define Am = 2 / (a^ + h'^^^'^rdr/ J aordr, where nam = / ^ cos{m(f))d(f) and irbm = 



- 12 - 



J T,sm{m(f))d(f) (cf. Boley et al. 2006 who use the volume density). In addition, define = 
Ylm=2 ^rn- As showu in Tablc 1. tends to increase with both resolution and decreased 
opacity. This is also seen in the Fourier component spectrum, shown in Figure 8, where 
squares indicate the 1024 simulations and crosses lines the 512 simulations. The 10^^ k 
simulation is not shown for readability, but if falls between the 10^^ and 10^"^ k, simulations 
and has the least amount of divergence between the 512 and 1024 simulations at large m. 
For all runs, the Fourier m — 1 through 6 dominate, and the power in the spectra 
drops off steeply around m ~ 10. The 1024 simulations have a shallower profile at large m, 
so resolution effects probably play a role for m > 10, even at resolutions of LMAX=512. 
The convergence of the Am profile for large m is a topic for future investigation (Steinman- 
Cameron et al., in preparation). 



4.2. Energy Budget 

The energetics of the disk simulations are shown in Figures 9 and 10, where cumulative 
energy loss by radiation (blue) and cumulative energy dissipation by shocks (red) are plotted. 
These quantities are integrated during the evolution, so they represent the actual heating 
and cooling in these disks. For both resolutions of the standard case, the shock heating 
clearly dominates over the cooling; the disk is heating up over the entire evolution. Energy 
is transported inefficiently for the highly optically thick disk, and the disk evolves like an 
adiabatic simulation. When the opacity is reduced by a factor of 100, the cooling becomes 
much more important than it is in the standard simulation, with radiative energy losses 
becoming comparable to heating by shocks. In addition, the shock heating rate becomes 
somewhat larger during the burst. The 1024 standard and 10~^ k, show additional heating. 
This heating is quite strong in the 1024 standard simulation, and it must be stopped even 
earlier than the 512 standard due to high temperatures, as discussed above. The 1024 10~^ 
K is very similar to its 512 counterpart. There is additional heating and cooling during the 
burst, but the heating and cooling curves for both resolutions roughly track each other. The 
adiabatic curve, which is only shown to 33 yr (1 orp), is difficult to see because it tracks the 
standard cases extremely closely. 

When the opacity is reduced by a factor of 10'^, the radiative cooling is even more 
efficient and surpasses shock heating. The 10^ opacity reduction shows the strongest shock 
heating and the fastest disk cooling. The 1024 10~^ k simulation curves track the 512 curves, 
with the burst contributing the most to the offset of the curves. Unlike the higher opacity 
simulations, for 10~^ k, there is less shock heating and less total coohng at higher resolution. 
The 1024 10~^ k, simulation roughly tracks the 512 counterpart, with the strongest deviation 
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near 66 yr (2 orp). The adiabatic curve is shown again as a reference. As one would expect 
from analytic arguments, the opacity has a profound effect on disk cooling and, consequently, 
on spiral shocks. 

The effect of opacity on energy losses is also demonstrated in Figures 11 and 12. In 
these figures, brightness temperature maps for the 1024 simulations are shown for similar 
times as in Figure 6, where = (ttJ+Zct)^/^ and /+ is vertical outward intensity as cal- 
culated by the ray solution in CHYMERA. In addition, we define a cooling temperature 
Tc = (/q°^ [—V ■ F] dz/aY^^, which is the temperature that corresponds to a given column's 
effective ffux if all of the energy were to leave the column vertically. If the column is being 
heated, is set to zero and appears as black in the images. 

As the opacity is lowered, the spiral structure becomes more clearly outlined, and the 
disk becomes brighter because the photons can leave from hotter regions. The brightness 
maps also demonstrate that sustained fast cooling by convection is absent. If hot gas near an 
optically thick midplane is quickly transported to altitudes where r ~ 1, convection can, in 
principle, enhance disk cooling. However, the efficacy of convective cooling is controlled by 
the rate at which that energy can be radiated from the photosphere of the disk. If locahzed 
convective flows were responsible for fast cooling in the optically thick disks, one would 
expect to find strong, localized Tb and T,. enhancements. For example, the standard and 
10~^ K, simulations would need to have regions as bright as the 10^^ k, simulation. There are 
thin regions along the spiral shocks in the 10^^ k simulation that show large T^, but they 
do not correspond to enhancements in T5 as well. This suggests that the energy from the 
spiral shocks is not transported efficiently to altitudes where it can be quickly lost from the 
disk. This assessment is supported by the net radiative heating in regions surrounding the 
strong shocks, e.g., the black outlines (areas experiencing net radiative heating) around the 
thin, hot spiral waves. The corresponding energy flux that is needed for convection to be 
sustained and to cool the disk does not occur (see also Boley et al. 2006, 2007b, Nelson 2006, 
& Rafikov 2007). 

The importance of the opacity for disk coohng is also shown in Figures 13 through 15. 
On these plots, three quantites are shown: the Toomre Q, the mass- weighted Fi, and the 
C = ^cooi^^//(ri) profiles. For the tcooi curves, the cooling times are calculated by dividing 
the azimuthally and vertically integrated internal energy by the azimuthally and vertically 
integrated radiative cooling rate for each annulus in the disk. The /(Fi) is the critical value 
of tcooi^j below which fragmentation is expected (Gammie 2001), for the corresponding Fi. 
Rice et al. (2005) demonstrated that /(7/5) i=s 12 and /(5/3) ~ 6. Based on these values, 
we assume /(Fi) pa — 23Fi -|-44 for the stabihty analysis presented here. One should keep in 
mind that /(Fi) is not a strict threshold and that the relation is approximate, especially for 
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simulations that permit the coohng time to evolve with the disk (Johnson & Gammie 2003; 
Clarke et al. 2008). Regardless, it serves as a general indicator for disk fragmentaiton. 

As Figures 13 through 15 demonstrate, the disk is approaching a state of constant Q 
for a wide range of radii, as expected in an asymptotic phase (Mejia et al. 2005). The 
mass-weighted average Qs between 3 and 6 AU are given in Table 1 and tend to decrease, as 
expected, when the opacity decreases. They do not seem to be greatly affected by resolution. 
The variation in the average Q with opacity is roughly consistent with the Aj^ measurements. 
One should keep in mind that the Q profiles represent snapshots of the disk, while the 
measurements arc time averages. The standard simulation is not shown in the stability plots, 
but the analysis of the 512 standard disk shows that it is highly stable against fragmentation. 

Because the cooling times fluctuate rapidly, we average the integrated cooling rates and 
the internal energies for the, roughly, 65- and 70/80-year snapshots. The cooling time profiles 
for the appropriate Fi only drop substantially below unity over extended regions where Q 
is high; these disks are stable against fragmentation. As the opacity is lowered, the coohng 
times decrease as well, with the 10~^ k simulation close to the fragmentation limit. It should 
also be noted that the only disk that behaves hke a constant Fi disk is the standard opacity 
simulation (not shown), and for all other simulations, a constant Fi is inappropriate (see 
Boley et al. 2007a). 

Mass in these disks is redistributed efficiently during the activation of GIs. Time- 
averaged mass fluxes are calculated by differencing the mass inside a cylinder at two different 
times. This represents the average mass flux as calculated by the second-order hydrodynam- 
ics scheme, and it is independent from the fluid element tracer. The inward and outward 
accretion rates vary, and in all simulations, can be well above 10~^ Mq yr~^ (Fig- 16). 

4.3. Shocks 

Due to our resolution, shocks are spread over scales of about 10^^ cm or larger, whereas 
ID calculations of chondrule-producing shocks give structures ~ lO^'^ cm (Desch & Connolly 
2002). To overcome this difficulty, we measure the pressure difference between the pre- and 
post-shock regions to calculate the Mach number M.. Using the pressure and temperature 
histories for each fluid element, a possible shock is identified whenever the dT/dt changes 
from negative to positive then back to negative. The pre-shock flow is taken to be the first 
sign switch, and the post-shock flow is the second sign switch. Let r] = (p2 — Pi) /pi be 
the fractional pressure change, where pi and p2 are the pre- and post-shock pressures, re- 
spectively, and where 7 is the average of the pre- and post-shock first adiabatic index. If 
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M"^ = (7 + 1)77/ (27) + 1 > 2, the event is counted as a shock. Once M. is determined, 
the pre-shock velocity Ui = Aicgi is calculated, where c^i is the pre-shock adiabatic sound 
speed of the gas. The shock strengths are also derived using the ratio of the post-shock to 
pre-shock temperatures: T2/T1 = - (7 - 1)) ((7 -1)M^ + 2)/ ((7 + ifM^). Both 

approaches should give the same answer whenever the shocks are cleanly defined and adia- 
batic. However, If radiative cooling becomes very efficient, then assuming adiabatic shock 
conditions may be incorrect. Furthermore, secondary waves or shocks (Boley & Durisen 
2006) may make identifying pre- and post-shcok regions very difficult. Neither method is 
clearly favored, and both will only provide crude estimates. 

We estimate a shock to be chondrule-forming when ui lies in a 1 km s^^ band between 
5 km s~^ and 11 km s~^ and when the pre-shock density is between logp( g cm~^) = -8.5 
and -10.5. This estimate is based on the Desch & Connolly (2002) ID calculations, similar 
to the relation used to locate chondrule-producing shocks in Figure 1. Chondrule-forming 
shocks should be thought of as having the potential to form chondrules, but may not yield 
chondritic material due to incompatible dust to gas ratios, too high or low cooling rates, and 
fractionation. 

Table 2 indicates shock information as extracted from the pressure histories of the fluid 
elements. Estimating shock strengths with the temperature ratio yields comparable numbers 
except for the chondrule-forming shocks (see below). Column two shows approximately the 
total number of detected shocks (TS) with Ai^ > 2, and column three displays the average 
number of shocks per fluid clement. Columns four (TS Mass) and flvc (CS Mass) show the 
total dust mass (see below) that goes through a shock with Ai^ > 2 and the total dust 
mass that encounters a chondrule-forming shock, respectively. To estimate how much dust 
is processed in shocks, we assign each fluid element a mass by calculating the total disk 
mass within some Ar and distributing that material evenly among all fluid elements in that 
Ar. A gas to solids ratio of 100 is assumed everywhere for the dust processing calculations. 
Although this is inconsistent with the growth and settling of solids assumed in six of the 
simulations, we do not worry about this detail inasmuch as the midplane will process more 
solids per shock and the high-altitude shocks will process less. Finally, column six shows the 
total number (Total FE) of fluid elements that remain after the 70 yr time period, i.e., the 
elements that were not accreted onto the star or fluxed into background density regions. 

We check whether the total number of detected shocks is reasonable by evaluating the 
expected number of shocks in a Keplerian disk with m-arm spiral waves. If the corotation 
radius for the pattern is r^, the number of shocks during some time period At for a fluid 
element orbiting at r is 



At. 



(6) 
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Evaluating equation (6) for 1000 fluid elements evenly distributed in annuli between 2 and 
8 AU, with = 4 AU, with m = 3, and with At = 60 yr yields approximately 8000 shocks. 
So the number of detected shocks in these simulations is reasonable. 

The total number of shocks per fluid element are roughly consistent for each simulation 
except for the 10~^ k runs. The more flocculent spiral structure in these simulations not 
only produces a larger number of shocks in the disk, but also creates candidate chondrule- 
producing shocks for four fluid elements; the 1024 10^'^ k, simulation has one fluid element 
that experiences a possible chondrule-forming shock. Several thousand of dust are 
pushed through shocks in each simulation, implying that almost all dust experiences a shock 
with > 2. In addition, a few of dust experience chondrule-forming shocks in the 
10~^ K runs. We remind the reader that these are only crude, order of magnitude estimates. 
When the temperature ratio is used, shocks with chondrule-forming strengths are undetected, 
but the total dust mass pushed through shocks remains the same to about one percent. 

5. Discussion and Conclusions 

In this section, we review the implications of the results of this study for disk frag- 
mentation, the effects of opacity on disk cooling, the FU Ori phenomenon, and chondrule 
formation. We remind the reader that the simulations presented here are meant to be a 
numerical experiment that explores the possible connection between chondrules, FU Ori 
outbursts, and bursts of Gl-activity. Moreover, this experiment provides a systematic study 
of the effects of opacity on disk cooling, which complements the Cai et al. (2006) metallicity 
study and provides a test bed for disk fragmentation criteria. 

5.1. Fragmentation 

After the onset of the burst, none of the disks fragments, and only the 512 10~^ k, 
simulation shows dense knot formation during the peak of the burst (10-30 yr). These knots 
do not break from the spiral wave even in the 1024 simulation, and so clump formation 
does not seem to be missed due to poor resolution. One should also keep in mind that for 
these simulations, the disk is flrst relaxed to equilibrium with standard opacity and then 
the opacity is suddenly dropped by a factor of 10*^ at t = 0. In effect, the dust settling is 
treated as instantaneous. Knot formation may not occur in a disk with more realistic settling 
timescales. As discussed below, this is also a caveat for the chondrule formation results. On 
the other hand, it does suggest that disk fragmentation might occur inside 10 AU under even 
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more extreme, but perhaps unphysical, conditions than we have modeled. 

The stabihty of these disks against fragmentation is supported by Figures 13 through 
15. The coohng rates for the standard and 10~^ k simulations are too low to cause disk 
fragmentation. The 10~^ and 10~^ k simulations do have areas where C, — icooi^//(ri) < 1, 
but the disks approach stability against GIs in those regions {Q > 1.7; Durisen et al. 2007). 
As noted in §4.2, ^ < 1 is not a strict instability criterion, but it does serve as an estimate 
for disk stability against fragmentation. For no simulation is C well below unity where Q is 
also < 1.7. 

Lowering the opacity increases the cooling rates. The 10~^ k simulations exhibit the 
most rapid cooling because the midplane optical depths are near unity, which results in 
the most efficient radiative cooling possible. Changes in dust opacity have a profound effect 
on disk cooling. Although not modeled, if the opacity were to continue to drop such that 
the midplane optical depth becomes well below unity, cooling would once again become 
inefficient. In addition, supercooling of the high optical depth disks (standard and 10~^ k) by 
convection does not occur. Based on our results, we believe that hydraulic jumps as a result 
of shock bores (Boley & Durisen 2006) rather than convection are a better explanation for 
the upwellings around spiral arms reported by Boss (2004a) and Mayer et al. (2007). These 
findings are consistent with analytic arguments by Rafikov (2005, 2007), with numerical 
studies of disk fragmentation criteria by Gammie (2001), Johnson & Gammie (2003), and 
Rice et al. (2005), and with global disk simulations where tests of the radiation algorithm 
and/or careful monitoring of radiative losses were performed (Nelson et al. 2000; Boley et 
al. 2006, 2007b; Stamatellos & Whitworth 2008). 

The radiative algorithm used in CHYMERA has passed a series of radiative transfer 
tests that are relevant for disk studies, including permitting convection when expected. 
Our conclusions regarding fragmentation are based on mulitple analyses: surface density 
rendering (Figures 4-7), an energy budget analysis (Figures 9 and 10), cooling temperature 
and brightness maps (Figures 11-12 ), and a cooling time stability analysis (Figures 13-15). 
It is also important to point out that these results do not contradict Stamatellos et al. (2007) 
or Krumholz et al. (2007), who find fragmentation in massive, extended disks (> 100 AU). 

It is also pertinent to demonstrate that the CHYMERA code can detect fragmentation 
when cooling rates are high and Q is low. Figure 17 shows a snapshot for a simulation similar 
to the 512 10~^ K simulation, but with the divergence of the fluxes artificially increased by 
a factor of two. In the normal simulation, kinks in the spiral waves form during the onset 
of the burst. As discussed above, the disk is very close to the fragmentation limit, but the 
knots do not break from the spiral wave. Figure 15 suggests, although for a later time, that 
increasing the cooling rates by a factor of two would drop C, well below unity in low-Q regions. 
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Figure 17 confirms that the disk fragments with such enhanced coohng. However, we remind 
the reader that the 10""^ k, simulation has the optimal optical depth for cooling (r ~ 1). 
We cannot imagine any physical process that could cause a factor of two enhancement in 
cooling. 

Three clumps form, one for each spiral wave, between 4 and 5 AU. The location of 
clump formation is consistent with the prediction by Durisen et al. (2008) that a spiral 
wave is most susceptible to fragmentation near corotation. One of the clumps survives for 
several orbits and eventually passes through the inner disk boundary. Resolution is always 
a concern for simulations. Because these results are consistent with analytic fragmentation 
limits and numerical fragmentation experiments, we conclude that CHYMERA can detect 
fragmentation at the resolutions employed for this study. 

As discussed above, when the opacity is abruptly decreased by a factor of 10^, the 
disk does approach fragmentation-like behavior. A similar effect is reported by Mayer et 
al. (2007), who find that their disk only fragments when the mean molecular weight is 
suddenly switched from /i — 2.4 to /i — 2.7. However, the simulation presented by Boley et 
al. (2006) was accidentally run with = 2.7, and Cai et al. (2006, 2008) purposefully ran their 
simulations at the high fi for comparison with the Boley et al. results. None of the studies 
reported disk fragmentation. This may indicate that when a disk approaches fragmentation 
shortly after a sudden switch in a numerical parameter, e.g., opacity here and fi in Mayer et 
al., the fragmentation may be numerically driven rather than physically, especially because 
sudden changes in cooling rates make disks more susceptible to fragmentation (Clarke et 
al. 2008). Regardless, the results do indicate that disk fragmentation by Gls may be possible 
under very extreme conditions. Whether such conditions are physically possible or realistic 
remains to be shown. Based on our results here and in earlier papers, disk fragmentation 
for r < 10s of AU appears to be a yet-unproven exception rather than the norm. 



5.2. FU Ori Outbursts 

In these simulations, the strong bursts of GI activity provide high mass fluxes (M > 
10~^ Mq yr~\ see Fig. 16) throughout each disk. Even though corotation is at r ~ 4 AU 
for the major spiral arms, the 2 AU region of the disk is strongly heated. Fluid elements 
approach peak temperatures of T ~ 1000 K in all simulations (Fig. 18). It is not difficult to 
speculate that if a larger extent of the disk were modeled, the temperatures due to shocks 
would be large enough to ionize alkalis thermally and possibly sublimate dust (1400-1700 
K). If such a condition is met, then an MRI could activate and rapidly carry mass into the 
innermost regions of the disk. From these simulations it appears to be plausible, but by 
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no means proven, that a burst of GI activity as far out as 4 AU can drive mass into the 
inner regions of the disk and create strong temperature fluctuations, which may then be 
responsible for a thermal instablility. Although we are simulating very massive disks, we 
note that such systems may exist during the Class 1 YSO phase. Additional studies need 
to be conducted, preferably with a self-consistent buildup of a dead zone, to address the 
efficiency of this mechanism in low mass disks. 

There are at least three observable signatures for this mechanism. First, if a Gl-bursting 
mass concentration at a few AU ultimately results in an FU Ori phenomenon, then one would 
expect to see an infrared precursor from the Gl burst, with a rise time of approximately tens 
of years. Second, one would also expect for a large abundance of molecular species, which 
would normally be frozen on dust grains, to be present in the gas phase during the infrared 
burst due to shock heating (Hartquist 2007, private communication). Third, approximately 
the first ten AU of the disk should have large mass fiows if the burst takes place near r ~ 4 
to 5 AU. We speculate based on these results that if the burst were to take place at 1 AU, 
then high outward mass fiuxes should be observable out to a few AU. 



5.3. Dust Processing 

Each simulation shows that a large fraction of material goes through shocks with A4^ > 
2. Although these shocks are weak, their abundant numbers may result in the processing 
of dust to some degree everywhere in the 2-10 AU region. In fact, such processing may be 
necessary for prepping chondritic precursors for strong-shock survival (Ciesla 2007, private 
communication) . 

Based on the arguments presented in §2, the intent was to produce spirals with large 
pitch angles by constructing a disk biased toward a strong, sudden Gl-activation near 5 AU. 
Even with this bias, the pitch angles of the spirals remain small, with i ^ 10°. Why are 
the pitch angles so small? According to the WKB approximation, coti =| krv/m \ (Binney 
& Tremaine 1987), where A:,, is the radial wavcnumber for some m-arm spiral. The most 
unstable wavelength for axisymmetric waves (Binney & Tremaine 1987) is roughly ~ 
2tiCs/Qk ^ 27rh/Q for disk scale height h. For a disk unstable to nonaxisymmetric modes, 
Xu corresponds to some m-arm spiral (e.g., Durisen et al. 2008). By relating kr — 27r/3m/A„, 

cotip^pQr/h (7) 

in the linear WKB hmit (| krv/m |» 1), where /3 is a factor of order unity that relates 
A„ to the m-arm spiral. Because (5Qr/h ~ 10 in gravitationally unstable disks, the linear 
WKB analysis may be marginally applicable. Equation (7) predicts that linear spiral waves 
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in these disks should have pitch angles i f=:i 6 to 11° for /3 between 1 and 1/2, respectively. It 
appears that this estimate for i extends accurately to the nonlinear regime, a result we did 
not expect. 

The GI spirals are efficient at heating the disk and transporting angular momentum, 
but not at producing chondrules. Only for the 10""^ k simulation are multiple candidate 
chondrule-producing shocks detected, and there is only one chondrule-producing shock in 
the 1024 10^^ K, simulation. These low opacity disks have relatively flocculent spiral mor- 
phologies, including kinks in spiral waves. The 10^"^ k, simulations arc on the verge of 
fragmentation. When interpreting these rcsTilts, it is important to remember that the opac- 
ities were lowered abruptly and in a quite unpliysical way. In addition, these detections are 
based on the fractional pressure change r]. If the temperature change is used, no chondrule- 
producing shocks are detected. This discrepancy may be due to efficient radiative cooling, 
which may make the assumption of adiabatic shock conditions incorrect, and/or to confusion 
with additional waves induced by shock bores. For the rest of this section, we assume that 
the pressure difference is the reliable shock identifier in order to discuss the implications of 
detecting chondrule-forming shocks in these disks. A thermal and spatial history for a fluid 
element that experiences a possible chondrule-forming shock is shown in Figure 19. 

To estimate the occurrence of a chondrule-forming shock, we employ a generous Ui- 
p criterion (Fig. 20). We do not take into account the optical depth criterion of Miura & 
Nakamoto (2006) on grounds that the large scale over which these shocks take place can allow 
for chondrules to equilibrate with their surroundings (Cuzzi & Alexander 2006). However, 
one should be aware that the optical depth criterion of Miura & Nakamoto will likely exclude 
all chondrule-producing shocks inside 4 AU in the Ui-p plane for these simulations. All 
candidate chondrule-forming shocks occur between r ~ 3 and 5 AU and at altitudes that are 
roughly less than a third of the gas scale height. Because a large degree of settling is assumed, 
these shocks are located in regions that may be consistent with the dusty environments in 
which chondrules formed (Wood 1963; Krot et al. 2005). 

We estimate that ~ 1 of dust is processed through chondrule-forming shocks. Be- 
cause these shocks are hmited to the onset of the GI burst, a fewxlO of dust would be 
processed in the 10~^ k disks if they went through about ten outbursts. If more outbursts 
take place than those that lead to an FU Orionis event, then more chondritic material could 
be produced. In order to produce these shocks, the disk was pushed toward fragmentation 
by suddenly dropping the opacity by a factor of 10^. We conclude that for bursts of GIs near 
4 or 5 AU to produce chondrules, the disk must be close to fragmentation. 
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6. The Unified Theory 

The hypothesis behind the work summarized here is that bursts of GI activity, dead 
zones, the FU Ori phenomenon, and chondrule formation are Unked. The general picture 
is that mass builds up in a dead zone due to layered accretion and that GIs erupt once Q 
becomes low. This activation of the instability causes a sudden rise in the mass accretion 
rate that heats up the disk inside 1 AU to temperatures that can sustain thermal ionization 
and an MRI. The MRI shortly thereafter activates the thermal instability (Bell & Lin. 1994; 
Armitage et al. 2001; Zhu et al. 2007). In addition, the strong shocks process dust and form 
chondritic material in the asteroid belt and at comet distances. 

Given the results of our simulations, we think the Boss & Durisen conjecture that dead 
zones bursting at r > 4-5 AU can produce chondrules is unlikely^ unless the disks are on the 
verge of fragmentation. Because the conditions for fragmentation inside r ~ 10s of AU are 
very difficulty to achieve, we conclude that this is not viable as a general chondrule-formation 
mechanism. As a result, we return to Figure 1. Based on our simple, analytic argument 
in equation (1), chondrule formation can take place in the asteroid belt and at cometary 
distances for Gl-bursting dead zones between about 1 and 3 AU, even with i 10°. As 
suggested by Wood (2005), chondritic parent bodies may be representative of temporal as 
well as spatial formation differences. Bursts that occur roughly between 1 and 3 AU seem to 
be able to accommodate this scenario, and could drive the FU Ori phenomenon (Armitage 
et al. 2001). As the disk evolves, the location of the dead zone can vary, with multiple bursts 
occurring at a few AU. Accretion rates between 10~^ and 10~^ Mq yr~^ could build up a 
0.01 Mq dead zone every 10'^ to 10^ years, respectively. Some of these bursts may drive 
the FU Ori phenomenon and produce chondrules, but others may only be responsible for 
chondrule-formation events. In this scenario, the disk need not fragment, and there should 
be between about 10 and 100 separate chondrule-formation events. Future work is required 
to ascertain the plausibility of this modified version of the unified theory. 
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Table 1: Simulation information. The name of the simulation indicates the fraction of the 
standard opacity that is used during the evolution of the disk. Qav is the average Q between 
3 and 6 AU for a snapshot near 70/80 yr, depending on the simulation (see text). The 
measurement is the sum of the time-averaged Fourier components between roughly 55 and 
80 yr. 



Sim. Name Resolution r, 0, z Duration Qav (^+) 



Adiabatic 


256, 512, 


64 


33 yr 






512 Standard 


256, 512, 


64 


80 yr 


1.7 


1.0 


1024 Standard 


512, 1024, 


128 


55 yr 






512 10-2 K 


256, 512, 


64 


110 yr 


1.7 


1.3 


1024 10-2 K 


512, 1024, 


128 


110 yr 


1.8 


1.6 


512 10-3 K 


256, 512, 


64 


110 yr 


1.5 


2.1 


1024 10-3 K 


512, 1024, 


128 


110 yr 


1.5 


2.2 


512 10-^ K 


256, 512, 


64 


110 yr 


1.4 


2.1 


1024 10-^ K 


512, 1024, 


128 


110 yr 


1.4 


2.4 



Table 2: Shock information for the time period between 10 and 70 yr. All estimates are based 
on the fractional pressure change r) (see text). TS indicates the total shocks encountered by 
all fluid elements. TS/FE gives the average number of shocks for a fluid element. TS Mass is 
a rough estimate of the total dust mass pushed through shocks in each disk, and CS Mass is 
the total dust mass that encounters a chondrule-forming shock. Finally, Total FE indicates 
the number of fluid elements that remain in the simulated disk at the end of the 60 yr time 
period. The 1024 standard simulation is omitted because it was stopped after 55 yr. 



Sim. 


TS 


TS 
FE 


TS Mass (Me) 


CS Mass (Me) 


Total FE 


512 standard 


8700 


9 


4(3) 





992 


512 10-2 K 


8600 


9 


4(3) 





990 


1024 10-2 K 


9500 


10 


4(3) 





962 


512 10-3 n 


8300 


9 


4(3) 





978 


1024 10-3 K 


9500 


10 


4(3) 


< 1 


964 


512 10-^ K 


9300 


10 


4(3) 


3 


939 


1024 10-^ K 


11000 


12 


5(3) 


2 


913 
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Fig. 1. — Expected shock speeds ui based on equation (1). In the legend, 
10) = Ui{rp = 1 AU, i = 10°) for corotation radius of the spiral pattern and pitch 
angle i. The colored regions highlight where chondrule formation is expected based on the 
Desch & Connolly (2002) shock calculations (see text), with the yellow region appropriate 
for a MMSN density distribution and the orange for the same density distribution but with 
10 X the mass. Shocks that occur inside corotaton for i = 10° will not produce chondrules 
between 1 and 5 AU. However, chondrules can be produced by these low pitch angle spirals 
in shocks outside corotation. If the pitch angle is fairly open, such as i ~ 30°, then a spiral 
wave with a corotation near 5 AU can produce chondrules in the asteroid belt and at comet 
distances. 
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Fig. 2. — Surface density profiles for the ICs with (sohd) and without (dashed) the density 
enhancement. The peak enhancement occurs near log(4.5 AU) ~ 0.65. 




Fig. 3. — Volume density contours for the initial model with (right) and without (left) the 
density enhancement. The normalization ro = 10 AU. 
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Fig. 4. — Surface density plots for the 512 simulations around 33 yr. 
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Fig. 5. — Surface density plots for the 512 simulations around 77 yr. 
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Fig. 6. — Surface density plots for the 1024 simulations around 33 yr. 
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Fig. 7. — Surface density plots for the 1024 simulations around 78 yr. The 1024 standard 
simulation is only run to about 55 yr (see text). 
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Fourier Component m 



Fig. 8. — Am profiles for the the 512 simulations (crosses) and the 1024 simulations (squares). 
The standard simulation does not have a 1024 counterpart because it was stopped at about 
55 yr (see text). These profiles demonstrate that the low-order structure dominates in the 
disk, with the Fourier components dropping off quickly near m ~ 10. The lower opacity 
simulations have stronger A^s, but there do seem to be some resolution effects present at 
large m. The 10"'^ k simulation is not shown for readability; it falls between the 10~^ and 
10~^ K, simulations and the 512 and 1024 resolutions follow each other closely. 
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Fig. 9. — Cumulative energy losses by radiation (blue) and heating by shocks (red) for the 
standard and 10~^ k simulations. The standard simulations are stopped early because the 
temperatures in the disk become too high to evolve due to inefficient cooling. The shock 
heating in the adiabatic (no cooling) simulation is indicated by the dark, dot-dash line. It is 
difficult to see because it follows the standard curves closely. For clarification, a A is shown 
adjacent to the cooling curves and a F is shown next to the heating curves. 
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Fig. 10. — Similar to Figure 9, but for the 10 ^ and 10 ^ k simulations. 
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Fig. 11. — Th and Tc maps for the 1024 standard and 10 ^ k simulations. As the opacity is 
lowered, the disks become much more efficient at cooling. See also Figure 12. 
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Fig. 12. — Same as Figure 11, but for the 1024 10 ^ and 10 ^ k simulations. 
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Fig. 13. — Stability assessment for the 1024 10~^ k simulations. The ordinate shows values 
for Q, Fi, and t cooi^/ fix i) as a function of r, where /(Fi) ^ — 23Fi + 44. litcooi^/ fO^i) > 1, 
then the disk is stable against fragmentation. 



- 41 - 




- 42 - 




Fig. 15. — Same as Figure 13, but for the 1024 10 ^ k simulation. 
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AU 

Fig. 16. — Time-averaged mass fluxes during the 0-55 yr period for tlie 512 Standard and 
tlie 1024 reduced opacity simulations. For each simulation, the mass fluxes are at FU Ori 
outburst levels. Positive values here correspond to net inflow, and negative values net outflow. 
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Fig. 17. — Similar to the 512 10 k simulation, but with the cooling artificially enhanced 
by a factor of two. The fragmentation is consistent with analytic predictions. 
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Fig. 18. — Initial r and minimum, average, and maximum temperature plots for fluid el- 
ements in the 512 Standard and the 1024 reduced opacity simulations. The temperature 
variations are quite large, approaching maximum temperatures of 1000 K. The values on the 
abscissa are due to the fluid elements that were lost from the simulated disk. 
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Fig. 19. — A thermal and spatial history for a fluid element in the the 1024 10~ k simu- 
lation. The fluid element shows strong variations in radius and height as it passes through 
shocks, with a net inflow due to gravitational torques. Strong variations in the pressure and 
temperature are also present. The strong shock near 30 yr is a potential chondrue- forming 
shock based on the pressure proflle. However, the temperature peak is not as high as one 
would expect for the same pressure jump in a ID adiabatic shock. This discrepancy may be 
due to efficient radiative cooling or to additional wave effects caused by shock bores. 
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Fig. 20. — Shocks on the ui-p plane for the 1024 10~^ k simulation, as determined from 
the fractional pressure change. The green area shows the chondrule-forming region, roughly 
based on the results of Desch & Connolly (2002). The green strip should not be thought 
of as definitive, and slight changes in the strip's location could identify more chondrule- 
forming shocks, particularly at low ui. If the pre- and post-shock temperature ratio is used 
to determine mi, no chondrule-forming shocks are detected and most shocks have a mi < 4 
km s~^. 



